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EXECUTIVE  SUMMARY 


Stress  analysis  of  composite  laminates  with  open  and  filled  holes  is  an  important  issue 
due  to  the  broad  range  of  composite  bolted  joining  applications.  The  inherent  difficulty  of 
the  three-dimensional  analysis  of  these  structures  is  due  to  singular  stress  behavior  near  the 
laminate  edge  and  ply  interfaces.  Numerical  stress  analysis  methods  such  as  commercial 
FEM  are  known  to  produce  mesh  sensitive  results  and  boundary  condition  errors  in  these 
regions.  This  report  reflects  further  development  of  the  combined  analytical  asymptotic 
solutions  and  B-spline-based  numerical  approximation  to  overcome  these  difficulties.  The 
method  of  superposition  of  a  hybrid  and  displacement  approximation,  developed  in  a 
previous  report  and  demonstrated  for  simple  laminates  with  open  holes,  is  extended  for 
singular  full-field  stress  analysis  in  multi-interface  practical  composites  under  mechanical, 
thermal  and  bearing  loading.  The  displacement  approximation  is  based  on  polynomial 
B-spline  functions.  This  method  provides  determination  of  the  coefficient  of  the  singular 
term  along  with  convergent  stress  components  including  the  singular  regions.  Reissner’s 
variational  principle  was  employed.  Uniaxial  loading  and  residual  stress  calculation  were 
considered  for  a  quasi-isotropic  IM7/5250  [45/90/-45/0]s  laminate.  A  convergence  study 
showed  that  accurate  values  of  the  coefficient  of  the  singular  term  of  the  asymptotic  stress 
expansion  could  be  obtained  with  coarse  out-of-plane  and  in-plane  subdivisions.  The 
interaction  between  singular  terms  on  neighboring  interfaces  was  found  to  be  important  for 
convergence  with  coarse  subdivisions.  Converged  transverse  interlaminar  stress  components 
as  a  function  of  distance  from  the  hole  edge  were  shown  for  all  examples.  The  hybrid 
approximation  developed  for  the  open-hole  problems  was  extended  to  rigid  fastener  hole 
problems.  A  near-singular  term  of  the  asymptotic  expansion,  present  in  the  filled-hole 
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problem  only,  was  taken  into  account  in  order  to  obtain  a  convergent  solution  with  a  coarse 
spline  approximation  mesh.  A  bearing-loaded  [45/-45]s  laminate  was  considered,  and  the 
multiplicative  factors  of  the  singular  stress  terms  at  the  hole  edge  and  ply  interface  were 
obtained . 
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1.  INTRODUCTION 


A  hybrid  approximation  model  was  proposed  in  [1,2]  for  singular  fiill-field  stress 
calculation  in  laminates  with  open  holes.  The  model  was  developed  for  single  interface 
laminates  under  mechanical  loading.  In  the  present  report  this  model  is  extended  for  stress 
analysis  in  practical  laminates  with  multiple  interfaces  under  thermal  and/or  mechanical 
loading.  Another  extension  reported  below  considers  rigid  fastener  hole  laminates  loaded  in 
bearings. 

The  model  developed  is  based  on  Reissner’s  variational  principle  and  is  intended  to 
reflect  the  singularities  which  arise  at  each  interface  at  the  boundary  of  the  hole.  Hybrid 
approximation  functions  to  be  developed  have  the  following  characteristics: 

(1)  They  include  the  asymptotic  solution,  thus  representing  the  directional 
nonuniqueness  of  the  solution.  It  is  only  in  this  manner  that  one  can  embed  the  proper 
singular  field.  The  fact  that  the  asymptotic  solution  results  from  the  three-dimensional 
problem  by  truncating  the  spatial  derivatives  in  the  circumferential  direction  [3]  will  be  used 
to  constract  hybrid  stress  functions. 

(2)  Two  independent  (B-spline)  displacement  functions  are  considered.  One  is  related 
to  the  regular  and  the  other  to  the  singular  portion  of  the  stress  field.  It  is  undesirable  to  use 
the  asymptotic  displacement  functions  in  the  displacement  approximation,  since  the 
calculation  of  their  derivatives  in  the  circumferential  direction,  required  in  the  variational 
formulation,  is  only  possible  numerically.  It  is  assumed  that  the  displacements  related  to  the 
singular  stresses  will  also  be  approximated  with  splines.  Thus  the  approximations  of  stresses 
and  displacements  are  made  independently. 
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2.  PROBLEM  STATEMENT 


Consider  a  rectangular  N-layer  laminate  built  of  orthotropic  layers  with  length  L  in 
the  x-direction,  width  A  in  the  y-direction,  and  thickness  H.  Individual  ply  thicknesses  are 
hp=z®-z"^‘’,  where  z=z^‘’^  and  z=z'‘^‘^  are  upper  and  lower  surfaces  of  the  p-th  ply,  respectively. 
The  origin  of  the  x,y,z  coordinate  system  is  in  the  lower  left  comer  of  the  plate,  as  shown  in 
Figure  1.  A  circular  opening  of  diameter  D  with  the  center  at  x=x^  and  y=y^  is  considered. 
Uniaxial  loading  is  applied  via  displacement  boundary  conditions  at  the  lateral  sides  (x=0,L): 

-  Wjc  (0,  y,  z)  =  Ux  (L,  y,  z)  =  Ur., 

(1) 

Uy(0,y,z)  =  Uy(L,y,z)  =  0, 

The  transverse  displacement  is  not  constrained  aside  from  a  rigid  body  constant,  and 
the  respective  traction  T2=0  is  prescribed.  The  edge  of  the  opening  is  part  of  the  traction 
prescribed  loading  boundary  S,.,  so  that 

Oijn.  =r,(x,.),x,GSj., 

or  mixed  boundary  conditions  in  the  case  of  the  rigid  fastener 

M,  =  0,  (JijTi  j  =  Ti  ix. ),  i  r,  Xi  e  Sj. 

where  tractions  7^.  are  known.  Indices  i,j=l,2,3  correspond  to  directions  x,y,z,  or  r,0,z, 
respectively.  The  lateral  edges  y=0,A  and  the  surfaces  z~0,H  also  belong  to  S,..  The 
constitutive  relations  in  each  ply  are  as  follows: 

^ij 

where  and  are  elastic  moduli  and  thermal  expansion  coefficients  of  the  p-th 

orthotropic  ply,  and  AT  is  the  temperature  change.  A  cylindrical  coordinate  system  r,  6,z  with 
the  origin  at  (jc^  y^,0)  is  introduced. 
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Figure  1.  Plate  with  the  Filied  Hole  and  Coordinate  Systems. 


x  =  rcosd  +  X  ,y  =  rsin0  +  y  ,z  =  z. 
c  c 

According  to  the  asymptotic  analysis  performed  in  earlier  papers  [3,4],  the  stresses  in 
the  vicinity  of  the  hole  edge  and  the  interface  between  the  p  and p+l  plies  are  thought  to  be 
of  the  form: 

Gy  =  ^  (6)7] (W^^)  +  bounded  function . 

A 


The  first  terms  represent  the  unbounded  (0<Re(A,)<l)  singular  stress  terms,  where  Ti  and  \|/ 
are  local  coordinates  introduced  in  the  cross  section  6=const  (Figure  l.b)  according  to 
equation 

T]  cosy/ =  r  -  D / 2,  T]  sinyr  =  z  -  z^’’^ .  (2) 

The  asymptotic  solution  is  normalized  similar  to  so  that 

K  (d)  =  lim[?7 (D / 2  +  7] cosi^, z , 0)] . 

The  singular  term 

is  a  solution  of  the  asymptotically  derived  two-dimensional  problem  and  satisfies  neither  the 
three-dimensional  equilibrium  nor  compatibility  equations  in  any  finite  volume  around  the 
curved  boundary.  However  it  does  satisfy  the  traction-free  boundary  conditions  on  the  hole 
edge  as  well  as  the  interface  continuity  conditions  in  the  limit  7]  0 .  We  shall  consider  one 

singular  term  at  each  angular  location  at  each  interface  at  the  open-hole  edge,  with  the 
extension  to  an  arbitrary  number  of  terms  at  the  same  singular  point  being  straightforward. 
The  displacement  components  and  bounded  portion  of  the  stresses  will  be  approximated  by 
using  cubic  spline  functions,  and  Reissner’s  variational  principle  will  be  applied  in  order  to 
obtain  Kp(0)  and  the  unknown  spline  approximation  coefficients.  The  asymptotic  solution 
near  the  orthotropic  ply  interface  and  the  hole  edge  will  be  considered  next. 
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3.  ASYMPTOTIC  SOLUTION 


Consider  a  region  around  the  hole  edge  and  the  interface  between  plies  p  and  p+1.  A 
local  coordinate  system  77,  y/is  introduced  in  the  radial  cross  section  0=const  according  to 

7t  7Z 

equation  (2).  In  this  coordinate  system  0<\f/<—  in  the  upper  ply  and - <\{/<0  in  the 

2  2 

lower  one.  For  an  arbitrary  function  F: 


where 


,  ^  dF  I  dF  .  ^  „  dF  .  1  dF 

AfF  = — cosy/ - smyr,  A^F  =  — smy/  + - cosi/^. 

dr]  7]dyf  dr]  t]  dy/ 


In  Cartesian  coordinates  the  derivatives  can  be  calculated  as 


^  =  (cos0)A^F- 
ax 

^  =  {smd)AfF+- 
dy 


sin0  dF 

DI2+T]cosy/  do  ’ 
cos  6  dF 

D/l+rjcosy/  dO  ’ 


(3) 


If  77  -»  0 ,  then  the  first  terms  in  the  right-hand  side  of  both  equations  (3)  are  of  the  order 
F/77,  and  the  second  terms  are  of  order  F.  Thus  for  small  77  the  expressions  for  derivatives  (3) 
are  truncated  by  retaining  only  the  first  terms,  so  that: 
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dx 

dy 

(dF\ 

dz 


dx 


dF 


=  — COS0  =  h.fF  cosd, 

je  *■  ' 


^dF^ 


[  dy 


Je 


dF 

=  — sin0  =  A*Fsin0,  77 ->0 
dr 


(4) 


)e 


=  ^  =  A„F 


where  the  notation  (  )g  implies  truncation  by  deletion  of  the  0  derivative.  Although  these 


equations  are  exact  in  the  limit  t]—^0  only,  we  shall  use  these  derivatives  through  a  finite 
volume  to  construct  the  hybrid  approximation.  Under  these  assumptions  the  Navier 
equations  of  a  given  ply  will  simplify  to: 


(AAA+BA„A,+CA„AJ 


=0, 


where  the  3  x  3  matrices  A3>C  are  given  in  the  Appendix  and  depend  upon  the  elastic 
moduli  and  0.  Note  that  the  thermal  expansion  coefficients  do  not  enter  these  equations  due 
to  the  assumption  of  uniform  temperature  change.  The  solution  of  the  Navier  equations  can 
be  found  in  the  form 

“i=X  vf(77^,v^,0)  +  U.'’(77,v^,0) 

A 

where  C//"  (77,  is  a  particular  solution  and  ^  \f  is  the  homogeneous 

A 

solution.  Nonhomogeneous  boundary  conditions  resulting  from  nonzero  prescribed  tractions 
or  thermal  mismatch  can  be  satisfied  with  a  piecewise  polynomial  particular  solution  of 
appropriate  order,  provided  that  the  prescribed  tractions  are  smooth  and  bounded  within  the 
ply  though  perhaps  discontinuous  at  the  interfaces.  The  thermal  expansion  strains  are  not 
present  for  the  homogeneous  solution  and  will  be  determined  in  a  particular  solution.  We 
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shall  be  interested  in  the  homogeneous  solutions  vf ,  which  have  essentially  nonpolynomial 

format,  and  the  respective  tractions  satisfying  homogeneous  boundary  conditions.  The 

homogeneous  solution  of  the  displacement  field  can  be  written  as  [3]: 

6 

vf  =n^5^7,ifjf(sinv^  +  /tf  cosvr)\ 

*=1 

where  the  superscript  p  refers  to  the  ply  number.  It  will  be  understood  that  coefficients 
and  vectors  { }  are  constants  for  each  ply;  therefore  the  superscript  is  subsequently 
omitted,  unless  needed  for  clarity.  The  stresses  from  the  homogeneous  solution  are 


where  the  tmncated  derivatives  are  calculated  by  using  expressions  (4).  The  expression  for 

the  stresses  may  be  written  in  the  form 

6 

+  l^k  COSVA)^"' , 

k=l 

Coefficients  c,*’^ ,  defined  in  [3],  depend  upon  elastic  moduli  of  the  p-th  ply.  It  was 

taken  into  account  that  the  following  relationship  applies  for  the  differential  operators 

6 

A,  vf  =  +  cosi/r)^-', 

k=\ 

6 

K  vf  =  ^V^'^J^Yk^kii^W  + cosxjf)^-^ , 

*=i 

which  leads  to  the  following  characteristic  equation  for  obtaining  : 

det[AiJ.l +  Bjj^  +C]  =  0, 

where  { }  are  eigenvectors  of  the  characteristic  matrix; 
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[Afi^  +  +C]  d ^2  =^- 

Ak^. 

The  power  %  and  coefficients  Tk  are  defined  from  the  boundary  conditions  of 
displacement  and  traction  continuity  at  the  interface  between  plies  p  and  p+1 

\|r=0;  vf(77,O,0)  =  v,r'(T?,O,0),  (nA 6)  =  erf ‘(77,0,0),  i  =  x,y,z, 

and  the  traction-free  open-hole  boundary  conditions: 

O'r^(r?,y,0)  =  O,  i  =  r,e,z, 

=  0,  i  =  r,e,z , 

or  the  rigid  fastener  frictionless  contact  conditions: 

wr‘(T7,^,0)  =  cTr‘(r7,^,0)  =  O,  i  =  e,z, 

u^{'n,~,d)=a^i7]~,e)=o,  i=e,z. 

Nontrivial  solutions  of  the  homogeneous  boundary  value  problem  exist  for  discrete 
values  of  X  only.  Coefficients  are  obtained  to  satisfy  the  interfacial  and  hole  edge 
boundary  conditions.  We  shall  be  interested  in  the  solutions  when  o<  Re(A)  <1.  These  terms 
provide  unbounded  stresses  which  dominate  the  solution  for  small  t|.  In  the  context  of  the 
present  work,  laminates  with  more  than  one  interface  will  be  considered,  and  the  singular 
asymptotic  terms  will  be  used  at  each  interface.  For  convenience,  we  will  introduce  an 
analytical  continuation  of  the  asymptotic  displacements  into  all  plies  of  the  laminate.  Thus 
we  extend  the  definition  of  the  displacement  vector  vf  and  stress  components  af  for  the 
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asymptotic  solution  associated  with  interface  between  plies  p  and  p+1  depending  upon  the 
properties  of  the  ply  in  which  \|/  is  located' 


ri^'^Ytdu^^sini}/ +  cosii/)^,\lf  >  0 
i=l 

6  ’ 
(sini/r  +  cos\iff,yf  <  0 

*=i 


(5) 


and 


afj=i 


(sini/r  +  cosyrr  \\j/>0,q  =  p  + 


fc=i 


(sinv^  +  )tif  cosyf)^~\\if  <0,q  = 


k=l 


(6) 


These  functions  are  defined  through  the  entire  laminate  thickness.  The  stresses  afj  are 


calculated  from  strains  defined  by  the  truncated  derivatives  of  displacements  vf  using  the 
elastic  moduli  of  that  ply  where  the  stresses  are  evaluated.  The  asymptotic  solution  obtained 
between  plies  p  and  p+1  will  provide  a  nonzero  stress  contribution  in  a  large  area 
surrounding  the  singular  point  due  to  the  weak  character  of  singularity  A,- 1—0.05  (references 
[1-4])  between  similar  orthotropic  layers  with  different  fiber  orientations  for  the  open  hole. 


*  Note  that  we  have  refrained  from  using  indices  on  the  variables  Tl  an  Y  since  they  always  emanate  from  the 
singular  point  at  z=z^\  Also  X  has  been  given  no  subscript  as  we  only  use  the  lowest  possible  value  at  z=z^\ 
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4.  VARIATIONAL  FORMULATION 


Displacement  components  are  represented  as 

U;  =u.  +u., 


(7) 


where  the  displacements  u.  satisfy  boundary  conditions  (1).  The  term  m/  is  associated  with 
singular  stress  components  and  the  second  term  u-  with  bounded  stress  components.  The 
stresses  will  be  assumed  as 


c  =<7*^*  +c^., 

y  y  y  ^ 


(8) 


The  stresses  and  displacements  u-  are  independently  assumed.  The  stresses  Cy 
and  displacements  u'  are  related  as  follows 


^ij  ~  ^ijkl  ^klj^T^)  » 


(9) 


where  and  a^i  are  elastic  moduli  and  thermal  expansion  coefficients  of  the  q-th 
orthotropic  ply,  AT  is  the  temperature  change,  and 


%;)  =  2 


3m  3m. 
- !-  + - 

dxj  dXj  j 


(10) 


Reissner’s  variational  principle  dR  =  0  is  employed,  where  the  functional  R  is 


given  by  equation 


R  =  III (-*(<7, , AT-)  +  )dV-\\T,  u,ds. 


(11) 


and 


*(<7, ,  AD  =  AT,  }=  }' 
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The  surface  tractions  Tj  in  equation  (1 1)  are  considered  in  a  generalized  sense,  meaning  that 
on  the  portion  of  the  surface  S,,  not  belonging  to  the  contact  surface  between  the  fastener  and 
the  hole,  they  are  prescribed.  On  the  contact  surface  these  tractions  are  unknown,  and  T; 
denotes  the  Lagrangian  multiplier  introduced  to  apply  the  boundary  condition  u  =0. 
Substituting  equations  (7)  and  (8)  into  (11),  one  obtains  after  algebraic  manipulations  and  use 


of  Betti’s  law; 


V 

+  ul^,,hT)dV -\\tM  +<)ds.  ’ 


where 


W(£,,Ar)  -alATXs,,  -a^AD 

The  goal  of  the  formulation  is  to  treat  problems  when  the  external  tractions  Tj  and/or  the 
interfacial  tractions  cannot  be  approximated  pointwise  by  using  the  same  shape  functions  as 
the  displacement  components  or  their  derivatives.  Let  the  functions  Xm  be  the  basis  functions 
for  approximation  of  the  displacements  ,  and  the  functions  Ym-  those  for  the 

displacements  u- .  Then  the  following  systems  of  equations  will  be  obtained  by  taking  the 
variations  with  respect  to  the  unknown  approximation  coefficients 

III  -  < ad + = J| T,x^ds 

V  Sj 

III  \ciuK„ -aiAT,  +  <jf},,jdV=\jT,Y,dS 

V  St 

Both  sets  of  basis  functions  Xm  and  Ym  in  the  present  context  are  of  the  same  piecewise 
polynomial  nature.  Without  restricting  the  generality,  one  might  assume  that  they  are  the 


13 


same.  Indeed,  one  might  always  redefine  the  systems  of  basis  functions  {Xm}  and  {Ym}  as 


{Xm}  u  { Ym} .  In  this  case  the  equations  can  be  rewritten  as: 

JJj  \PluK.n  +«(U,  - A7-a«)]x.jdy  = 

V  St  (13a) 

III  =  JJJ  (T^X„,dV  (13b) 

V  V 

Integrating  by  parts  and  adding  the  equations  to  each  other,  one  obtains 

JJJ  -Ara')  +  CTf  )n,]x.<B=0 

V  Sj 

The  first  term  of  the  above  equation  contains  the  weak  form  of  the  equilibrium  equations  and 
the  second  term  -  the  weak  form  of  the  boundary  conditions.  They  are  interconnected, 
meaning  that  if  the  boundary  conditions  are  not  satisfied  in  the  weak  form,^  then  the  error  in 
the  equilibrium  equations  will  not  vanish  even  in  the  weak  form,  i.e.  it  will  not  be  orthogonal 
to  each  of  the  displacement  approximation  basis  functions.  In  the  next  section  a  form  of 

will  be  proposed  so  that  the  boundary  conditions  on  the  hole  edge  and  the  interfacial 
continuity  conditions  will  be  satisfied.  In  this  case,  provided  that  equations  (13a,b)  are  also 
satisfied,  the  stresses  (8)  will  satisfy  the  equilibrium  equations  in  the  weak  sense,  even  in  the 
vicinity  of  the  singularities. 


2 


The  error  is  orthogonal  to  each  of  the  displacement  approximation  basis  functions. 
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5.  HYBRID  APPROXIMATION 


Consider  the  exact  stresses  associated  with  the  displacement  wf  in  the  q-th  ply: 


The  thermal  expansion  term  is  not  included  with  the  singular  displacement  portion,  since  it 
was  accounted  for  in  (9).  The  stresses  resulting  from  displacement  field  m/  are  modified  to 
include  the  singular  asymptotic  stress  field  (6).  We  shall  calculate  the  strain  field  generated 
by  the  truncated  derivatives  of  the  displacements,  as  follows 


(dul  ^ 

+ 

ra«;  ^ 

0 

dx, 

<  '  / 

'e 


where 


dXj 


are  calculated  according  to  truncated  expressions  (4).  The  contribution  of  the 
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stresses  generated  by  these  strains  is  given  by  the  asymptotic  stress  field  (6),  so  that  the 
hybrid  stresses  associated  with  displacements  u’  are 


N-l 


/’=! 


where 


^ij  '  (14) 

The  unknown  functions  of  the  hybrid  approximation  (7)  and  (8)  are  the  displacements  u\ ,  wf 
and  coefficients  K^.  Practically  speaking,  we  expect  the  hybrid  approximation  to  be  needed 
only  in  the  vicinity  of  the  hole  edge.  Let  the  volume  F  be  bounded  by  the  hole  edge,  the  top 
and  bottom  surfaces  and  r=ro:  D/2<r„<min(L,A).  Inside  this  region  stresses  are  given  by 
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equations  (8),  and  outside  this  region  m/  =  0  and  cr,^.  =  Gy .  The  total  displacement  must  be 

continuous  at  the  internal  boundary  r=ro.  It  can  be  satisfied  easily  if  the  unknown 
displacement  approximation  functions  are  the  total  displacement  u.  and  u'  instead  of  u'  and 
m/  .  The  stress  approximation  (8)  can  be  rewritten  accordingly: 


where 


JV-l 


p=\ 


(15) 


(16) 
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6.  GOVERNING  EQUATIONS 

Taking  into  account  equations  (15)  and  the  change  to  independent  displacement 
functions  u.  and  m/  ,  equations  (13)  will  attain  the  following  form 

III  [q«(«,..,>-A7-as)k.j‘n'=Jj7;x.<is  (17) 

V  Sr 

JJI  [q««i,o,  =  III  XX,  (e)a|X„^,dV  (18) 

r  r  p=i 

where  are  based  on  truncated  derivatives  (4).  Equation  (17)  follows  from  (13a)  after 

substituting  (7)  and  allows  one  to  calculate  the  total  displacement  in  an  independent  problem 
under  the  given  traction  and  displacement  boundary  conditions  (1).  Several  steps  have  been 
undertaken  to  reduce  (13b)  to  (18).  It  has  been  considered  that  the  hybrid  approximation  is 
defined  inside  T  only  and  that  tmncated  derivatives  (4)  were  used  in  the  hybrid  stress 
approximation.  The  spatial  derivatives  appearing  in  (18)  are  also  truncated  derivatives, 
meaning  that  auxiliary  problem  (18)  is  a  two-dimensional  problem  with  the  0-coordinate  as  a 
parameter.  A  substitution 

N-\ 

P=1 

provides  a  solution  of  equation  (18)  for  arbitrary  Kp  if 

111  [q«»S),]x.j|.d''=|||  (19) 

r  r 

The  displacement  boundary  conditions  on  the  singular  components  must  exclude  rigid 
body  motion  and  be  consistent  with  (5),  i.e.  no  surface  displacements  can  be  prescribed  for 
on  Sj-  n  dr.  Equation  (19)  provides  the  weak  form  of  the  equality  of  tractions  aF.rij 

and  on  the  portion  of  dr  inside  the  laminate  (r=ro),  provided  that  no  surface 
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displacements  are  prescribed  for  there.  If  the  distance  between  this  boundary  and  the 

hole  edge  is  chosen  sufficiently  large,  then  the  resulting  traction  discontinuity  at  this 
boundary  for  the  stress  field  (15)  may  be  reduced  arbitrarily  by  increasing  the  numerical 
subdivision  for  the  auxiliary  problem.  The  following  boundary  conditions,  consistent  with 
(5),  were  imposed 

(Z) / 2, 0)  =  ul'” (D ,e)  =  0,  ul'" (r„z^'’\e)  =  0,  (20) 


The  stress  approximation  after  the  substitution  can  be  rewritten  as 


<T  =<T“ 

ij  ij 


AT-l 

P=1 


where 


5.? 


s,p 


(21) 
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7.  DETERMINATION  OF  Kp(0) 

The  error  in  boundary  conditions  near  the  singularity  from  (17)  is  a  result  of 
approximating  the  directionally  nonunique  singular  stress  field  by  polynomials.  The 
directional  nonuniqueness  means  that  for  j]  0,  the  stresses  aF.  may  tend  to  plus  or  minus 

infinity  depending  upon  \|/.  The  polynomial  approximation  provides  unique  and  finite  stress 
values  at  every  point  of  a  given  ply.  The  values  of  the  stress  intensity  factors  are  obtained  to 

■  N-\ 

enforce  single-valuedness  of  the  nonsingular  portion  t7,“  -  {d)s[j  (O)iafj  -  sfj )  at 

9=1, 

q*p 

the  singular  point.  Consider  the  interface  z=z''’*  between  the  p  and  the  p+i-st  ply.  Let 
be  the  interlaminar  traction  on  the  interface  the  vicinity  of  the  edge 

( (0)=  ±00).  Then,  at  the  hole  edge  six  traction  boundary  conditions  have  to  be 
simultaneously  satisfied 

Ti=Cijn';,Pi^'’\ri)=ayn^, 

Thus  one  stress  component,  namely  o„,  will  appear  in  two  equations  which  are 

These  equations  can  only  be  satisfied  if 


-  Kp  (77.-,  0)1  =  -  Kp  (e)a^  (77.0, 9)]  = 


These  equations  are  necessary  conditions  to  make  the  polynomial  part  of  the  stress  tensor 
single-valued.  It  is  worth  noting  that  if  instead  of  a  comer  one  has  a  crack,  i.e.  nf  =  -nf , 
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then  we  shall  have  three  pairs  of  single-valuedness  conditions,  one  per  each  stress  component 
in  the  plane  normal  to  the  crack.  It  will  require  three  stress  intensity  factors:  Mode  I,  n  and 
in.  In  the  case  of  a  hole  edge,  we  have  only  one  interlaminar  stress  component  which 
requires  the  single-valuedness  condition:  .  The  criterion  for  determining  KJ[ $)  will  be 

the  continuity  of  the  nonsingular  part  of  cr^  stress  component  at  the  singular  points,  i.e. 


=  K  (e)A 

L  y  2  j 

p 

_  y  2 

N-l 


+  K^(d)A 


9=1. 

g*p 


,  p  =  l,...,N-l 


(23) 


where  a[  ]  denotes  the  difference  of  the  bracketed  function  between  z^’+O  and  z®-0.  The 
nondiagonal  terms  in  the  right  side  of  equation  (23)  represent  the  influence  of  the 
singularities  of  adjacent  plies.  These  terms  will  become  small  if  the  subdivision  of  the  p  and 
p+1  is  dense,  since  the  nondiagonal  terms  contain  the  difference  between  the  asymptotic 
solution  and  its  polynomial  approximation  away  from  the  singularity.  However,  for  coarse 
subdivisions  the  contribution  of  the  adjacent  plies  is  significant  for  convergence. 
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8.  SPLINE  APPROXIMATION  OF  DISPLACEMENT  COMPONENTS 


The  X,  y  and  z  displacement  components  are  approximated  by  using  cubic  spline 
functions  in  curvilinear  coordinates.  The  total  displacement  is  approximated  as: 


M.  =CiXU’'  +(5.iM„XE^, 


(24) 


where  X  is  a  vector  of  three-dimensional  spline  approximation  basis  functions,  and  U;  are  the 
unknown  spline  approximation  coefficients.  The  nonsquare  matrices  Q  and  constant  vector 
E  are  defined  so  that  the  approximation  (24)  is  kinematically  admissible,  i.e.  it  satisfies 
boundary  conditions  (1)  for  arbitrary  coefficients  Uj.  Bold  type  will  be  used  to  distinguish 
vectors  and  matrices;  superscript  ^  means  the  transpose  operation. 

A  detailed  description  of  the  spline  approximation  procedure  and  the  properties  of 
spline  functions  are  given  by  larve  [3].  The  three-dimensional  spline  approximation 
functions  are  defined  in  curvilinear  coordinates.  The  x,y  plane  was  mapped  into  a  region 
defined  by  p,  (j) ,  where  0  < p  <  1  and  0  <(j) ^2k.  The  transformation  was  defined  as  follows: 


x  =  ^I\  (P)cos^  +  L  ■  i^(p)a(0)  + 
y  =  ~  Fj(p)sin  (j>+A-  F2(P)A^)  +  Jc 


(25) 


Functions  Fj  and  F2  were  defined  as 


Flip)  = 


l+KT-p.p  <  p^ 

(l  +  g  P/,)(l-p),Pfc  <  P  <  1 


F2ip)  = 


i-Ph 


0,p<p^ 

P-Ph  ,Pf,<P<l 

^-Ph 


Coordinate  line  p=0  describes  the  contour  of  the  hole,  and  the  coordinate  line  p=l  describes 
the  rectangular  contour  of  the  plate.  Parameter  k  defines  the  size  of  a  near-hole  region 
D/2  <r  <(1+k)D/2,  which  corresponds  to  0<p  <ph,  where  a  simple  relationship  between  the 
polar  coordinates  and  the  curvilinear  coordinates  p,  ^  exists: 
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(26) 


D  Dk  j  Cl  a 

r--  =  ~p  and  e=(t). 

The  width  of  this  region  is  typically  one  hole  radius,  i.e.,  =  1 .  Beyond  this  region  a 

transition  between  the  circular  contour  of  the  opening  and  the  rectangular  contour  of  the  plate 
is  performed.  Functions  a(<t))  and  |3(<t))  describing  the  rectangular  contour  of  the  plate 
boundary  were  given  by  larve  [3].  These  functions  are  introduced  so  that  parametric 
equations  x=a((l))+Xc,  y=P(<l>)+yc  describe  the  rectangular  contour  of  the  plate,  and  0  <{j)<(})^^^ 
corresponds  to  0<x<L,y=A;  (j)^^^<(|K(t)^^^  corresponds  to  x=0, 0<y<A; 
corresponds  to  0  <x<L,  y=0;  and  <(tK27C  corresponds  to  x=L,  0  <y<A. 

Sets  of  one-dimensional  cubic  spline  functions  were  defined  in  each  coordinate 
direction.  Subdivisions  were  introduced  through-the-thickness  of  each  ply  such  that  the  p-th 

ply,  which  occupies  a  region  <  z  <  ,  is  subdivided  into  np  sublayers.  A  set  of 

=  S/tp  +2N  +  l  basic  B-type  cubic  spline  functions  •^j(z)|^j  with  variable  defect  was 

built  along  the  z-coordinate  according  to  a  recurrent  procedure  given  by  larve  [3].  These 
cubic  splines  are  twice  continuously  differentiable  at  all  nodal  points  inside  each  ply  and  have 
discontinuous  first  derivatives  (the  function  itself  is  continuous)  at  the  ply  interfaces  to 
account  for  interfacial  strain  discontinuities.  This  is  a  necessary  condition  for  traction 
continuity  at  the  ply  interfaces.  Nodal  points  are  also  introduced  in  the  p  and  <|)  directions  as 
follows:  0  =  Pq  <  Pj  <  ....  <  p;„  =  1 ,  0  =  < ....  <(pi^  =2k.  The  subdivision  of  the  p 

coordinate  is  nonuniform.  The  interval  size  increases  in  geometric  progression  beginning  at 
the  hole  edge.  The  region  0  <p  <ph  in  which  the  curvilinear  transformation  is  cylindrical  is 
subdivided  into  m^  intervals  (mi<m),  so  that  p^=Pm  .  Sets  of  basic  cubic  spline  functions 
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r  r  •y:+3 

^  ,  -p  ^  along  each  coordinate  are  built  so  that  twice  continuous  derivatives  in 


each  node  are  provided.  Splines  along  the  (j)  are  periodic  at  the  ends  of  the  interval.  The 
vector  of  the  three-dimensional  spline  approximation  functions  was  defined  as  the  tensor 
product  of  three  one-dimensional  sets  of  splines: 

g  =  l  +  (j-l)N  +(i-l)N  k,  l  =  ,j  =  l,...,k,i  =  l,...,m  +  3. 

z  z  z 

The  components  of  vector  E  are  equal  to  1  or  -1  for  a  component  of  X  that  is  nonzero 
at  p=l,  <l)^^^<(tx(t)®  (x=0)  and  p=l,  <^<2k  (x=L),  respectively.  All  other  components  of 
the  vector  E  are  equal  to  zero.  The  boundary  matrices  are  obtained  by  deleting  a  number  of 
rows  from  the  unit  matrix.  The  deleted  rows  have  nonzero  scalar  product  with  E. 

The  region  T  of  the  hybrid  approximation  superposition  is  inside  the  region  in  which 
the  transformation  (25)  coincides  with  (26).  The  boundary  r=ro  is  defined  to  coincide  with 


the  radial  coordinate  line  so  that  r^=^(l  +  Kp^^),  where  <  mj .  A  reduced  set  of  splines 


in  the  p  -direction  (P)T=i^  defined  only  over  the  first  rrij  intervals  for  the  purpose  of 
efficient  solution  of  the  systems  of  equations  for  determining  .  It  is  built  exactly  the 

same  way  as  the  one  over  the  entire  interval.  It  can  be  shown  that  the  reduced  set  is  a  subset 
of  the  approximation  (24). 

The  approximation  of  the  singular  displacements  can  be  written  as 


u;-”  =CfX'Uf 


(27) 


where  Uf  (6)  are  unknown  coefficients  and 
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=  Rhp)Zj(z), 

q  =  l  +  (i-l)N  , 
z 

I  =  ,j  =  +3. 

Z  1 

Matrices  Cf  are  defined  to  satisfy  boundary  conditions  (20).  The  truncated  in-plane 


derivatives  in  equation  (19)  are  calculated  as: 


dy 


2  d 
= - sin^ — , 


Je 


Dk 


dp 
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9.  NUMERICAL  RESULTS 
9.1  Uniaxial  Tension  of  [45/907-45/0]^  Laminate 

A  [45/90/-45/0]j  quasi-isotropic  IM7/5250-2  laminate  was  considered,  where  the 
stacking  order  is  from  the  top  to  the  central  plane,  reading  from  left  to  right,  respectively. 
The  elastic  properties  of  the  unidirectional  ply  were  E,  =151  GPa,  £^=£3=9 .45  GPa, 
Gij=Gi3=5.9  GPa,  G23=3.26  GPa,  v,2=v, 3=0.32,  and  V23=0.45.  The  in-plane  dimensions  of  the 
plate  were  L=290  mm,  A=76  mm,  x^=L/2,  y=A/2,  D=12.5  mm  and  the  ply  thickness  h=1.34 
mm. 

The  average  applied  stress  was  calculated  as 

^0  =-^j^jy„(L,y,z)dydz  (28) 

Figure  2  shows  the  power  of  singularity  calculated  for  each  interface  as  a  function  of 
the  polar  angle;  it  varies  for  stresses  from  -0.01  to  -0.077  depending  on  the  angle.  Two 
mesh  densities  were  used  for  obtaining  the  coefficient  of  the  singular  term  under  the  uniaxial 
loading  boundary  conditions  (1).  Coarse  subdivision  was  defined  as  follows, 
z-coordinate:  1  sublayer  per  ply 

p-coordinate:  m=12, mo=8,Kp^  =  l,  q=1.2. 

6-coordinate:  48  equal  intervals. 

The  fine  subdivision  consisted  of: 

z-coordinate:  3  nonuniform  (1:2:1)  sublayers  per  ply 
p-coordinate:  m=24,  mg=20,  Kp^  =  l,  q=1.2. 

0-coordinate:  48  equal  intervals. 
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— — 0/-45  interface 
■  -45/90  interface 

90/45  interface 


Figure  2.  Power  of  Singularity  versus  Cylindrical  Angle  at  all  Interfaces  of  the  [45^/-45/0]s 
Laminate. 

The  coefficient  of  the  singular  stress  term  calculated  by  using  these  subdivisions  is  shown  in 
Figure  3  at  each  interface.  The  values  obtained  by  using  equation  (23)  and  the  values 
calculated  by  retaining  only  the  diagonal  terms  on  the  right-hand  side  of  equation  (23)  are 
compared.  For  the  coarse  subdivision  the  nondiagonal  terms  are  very  significant,  so  that  the 
values  obtained  by  neglecting  them  may  even  be  of  different  sign.  For  the  fine  subdivision 
the  two  results  are  almost  identical.  Indeed,  the  magnitudes  of  the  nondiagonal  terms  are 
determined  by  the  accuracy  of  approximation  of  the  singular  stress  term  by  the  polynomial 
approximation  one  or  more  ply  thicknesses  away  from  the  singularity.  The  density  of  the 
subdivision  through  the  ply  thickness  defines  this  accuracy.  The  values  of  the  K^,  Kj  and 
obtained  with  the  two  subdivisions  and  taking  into  account  the  nondiagonal  terms  are  close 
together.  Figure  4  shows  the  characteristics  of  the  singular  behavior  of  the  interlaminar  shear 
stresses,  namely: 
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--0--12  int.,  IsbI.  -  eqn.(23) 
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Cylindrical  angleO 

Figure  3.  Singular  Term  Coefficients  at  Different  Interfaces  for  the  [45/90/-45/0]s  Laminate 
under  Uniaxial  Loading. 
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Figure  4.  Transverse  Shear  Stress  Normalized  Singular  Term  Coefficients  at  Different 
Interfaces  for  the  [45/90/-45/0]s  Laminate  under  Uniaxial  Loading. 
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(Z)/2+77,z^^>,0)  =  is:  (1A6>) 

limr? {DI2+r],z^^\e)  =  Kal (1,0, 0)  ’ 

where  were  determined  by  the  fine  mesh  solution. 

Interlaminar  stress  components  at  a  0=60.3®  cross  section,  normalized  by  the  average  applied 
stress  (28),  will  be  considered  in  Figures  5-8.  The  transverse  radial  shear  component  on  all 
interfaces  calculated  by  using  the  fine  subdivision  is  shown  in  Figure  5.  Stresses  cy“rz  are 
displayed  in  Figure  5a.  At  the  0/-45  and  -45/90  interfaces,  we  observe  the  typical 
discontinuity  of  a“rz,  near  the  singularity,  whereas  the  interface  90/45  shows  relatively 
continuous  traction  up  to  the  singular  point.  This  is  in  agreement  with  Figure  4c,  showing 
very  small  amplitude  of  the  singular  term  at  the  latter  interface  compared  to  the  two  others  at 
0=60.3®.  The  stresses  calculated  according  to  hybrid  approximation  (21)  are  shown  in  Figure 
5b.  It  should  be  noted  that  the  closest  point  to  the  singularity  shown  in  this  figure  is 
r=D/2+0.002H.  The  traction  discontinuity  is  practically  indistinguishable  within  the  graphic 
resolution. 

The  transverse  normal  stresses  are  examined  in  Figures  6  and  7.  The  a“zz  stresses  on 
all  three  interfaces  are  shown  in  Figures  6a  and  6b.  The  difference  between  the  results 
obtained  with  two  subdivisions  is  clearly  observed  for  (r-D/2)/H<0.6.  The  values  obtained  by 
using  the  hybrid  approximation  are  displayed  in  Figure  7.  In  this  case  the  two  subdivisions 
give  practically  identical  results.  We  observe  (Figure  8)  a  small  absolute  value  of  K  (60.3®) 
at  the  90/45  interface.  This  led  to  an  inconclusive  Ozz  trend  definition  based  on  equation  (17) 
as  shown  in  Figure  6a.  The  hybrid  approximation  in  Figure  7a  shows  that  this  interface  is  in 
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Figure  6.  Transverse  interlaminar  Normal  Stress  at  6=60.3°  Obtained  by  Using  Displacement 
Approximation. 
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(b) 


Figure  7.  Transverse  Interlaminar  Normal  Stress  at  0=60.3°  Obtained  by  Using  Hybrid 
Approximation. 


0/-45  interface 


Figure  8.  Transverse  Interlaminar  Shear  Stress  at  6=60.3°. 
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compression  along  with  the  0/-45  interface.  The  -45/90  interface  exhibits  a  tensile  peel  stress 
singularity,  which  is  in  agreement  with  Figure  8. 

The  obtained  by  using  the  hybrid  approximation  with  the  two  subdivisions  are 

shown  in  Figure  8  for  completeness.  The  fine  and  the  coarse  subdivision  provide  very  close 
agreement  in  the  stress  values  which  is  indicative  of  convergence. 

9.2  Thermal  Stresses  in  a  [45/90/*45/0],  Laminate 

The  same  laminate  is  considered  under  a  uniform  temperature  change  of  AT=-167“C. 
This  temperature  drop  is  used  to  approximate  the  residual  stresses  generated  during  the 
processing  cool-down  phase.  The  displacement  boundary  conditions  (1)  are  modified  so  that 
the  edge  x=L  is  released  (zero  tractions).  The  coefficients  K^,  Kj  and  obtained  with  the 
coarse  and  fine  subdivisions  described  in  the  previous  section  are  shown  in  Figure  9.  Some 
differences  between  the  coarse  and  fine  mesh  analysis  results  can  be  seen  for  low  values, 
while  all  the  maximum  values  are  practically  identical.  The  values  of  the  coefficients  are 
positive  at  all  interfaces  in  all  circumferential  directions.  Therefore,  a  tensile  peel  stress  is 
expected  near  singularities.  Figure  10a  shows  the  transverse  normal  stresses  c"  obtained 
with  the  coarse  subdivision  based  on  displacement  approximation  in  the  cross  section 
0=157.8”,  where  the  maximum  value  of  K  occurs.  The  stresses  are  displayed  at  two 
interfaces  in  the  top  (solid  line)  and  bottom  (dashed  line)  plies.  The  normal  stresses  are 
seemingly  discontinuous  a  distance  of  approximately  0.6H  from  the  hole  edge.  Figure  10b 
shows  the  same  stresses  at  all  interfaces  calculated  using  the  hybrid  approximation  (21).  The 
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Figure  9.  Singular  Term  Coefficients  at  Different  Interfaces  for  the  [45/90/-45/0]s  Laminate 
under  Thermal  Loading. 


while  all  the  maximum  values  are  practically  identical.  The  values  of  the  coefficients  are 
positive  at  all  interfaces  in  all  circumferential  directions.  Therefore,  a  tensile  peel  stress  is 
expected  near  singularities.  Figure  10a  shows  the  transverse  normal  stresses  (7^  obtained 
with  the  coarse  subdivision  based  on  displacement  approximation  in  the  cross  section 
6=157.8°,  where  the  maximum  value  of  K  occurs.  The  stresses  are  displayed  at  two 
interfaces  in  the  top  (solid  line)  and  bottom  (dashed  line)  plies.  The  normal  stresses  are 
seemingly  discontinuous  a  distance  of  approximately  0.6H  from  the  hole  edge.  Figure  10b 
shows  the  same  stresses  at  all  interfaces  calculated  using  the  hybrid  approximation  (21).  The 
discontinuities  are  reduced  significantly.  It  is  interesting  to  point  out  that  the  <7^  stress  at  the 
90/45  interface  obtained  with  coarse  subdivision  in  Figure  10a  may  appear  to  have  a 
tendency  towards  -oo  as  one  approaches  the  interface.  The  hybrid  approximation  shows  a 
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Figure  10.  Transverse  Interlaminar  Normal  Stress  at  0=157.3°  Obtained  by  Using 
Displacement  Approximation  with  Coarse  Subdivision. 


sharp  change  of  sign  very  close  to  the  singularity.  This  trend  is  picked  up  by  the 
approximation  using  fine  subdivision  as  shown  in  Figure  1  la.  In  this  case  the  tractions  at  the 
bottom  and  top  surfaces  are  continuous  starting  from  0.2H  from  the  hole  edge.  The  hybrid 
approximation  (Figure  1  lb)  provides  traction  values  with  imperceptible  discontinuities 
everywhere.  The  hybrid  stresses  obtained  with  the  two  subdivisions,  Figures  10b  and  1  lb, 
are  practically  the  same. 

9.3  Bearing  Loading  of  a  [45/-45],  Laminate 

A  square  [45/-45]s  plate  is  considered.  The  geometric  properties  are  L=A=508  mm, 
D=50.8  nun,  ply  thickness  h=5.1  mm,  Ei=138  GPa,  E2=E3=14.5  GPa,  Gi2=Gi3=G23=5.9  GPa 
and  Vi2=Vi3=V23=0.21.  The  displacement  boundary  conditions  were  applied  so  that 
uo/L=0.001,  The  average  applied  stress  ctq  was  calculated  afterwards  according  to  equation 

(28)  and  used  to  normalize  the  results.  Formulation  and  solution  of  the  three-dimensional 
contact  interaction  problem  based  on  B-spline  displacement  approximation  is  given  in 
reference  [5],  Two  subdivisions  were  used  for  the  convergence  study.  The  0  coordinate  in  all 
cases  was  uniformly  divided  into  48  intervals.  The  “coarse”  subdivision  consisted  of  ns=l- 
one  sublayer  per  ply  in  the  z-direction,  and  a  total  of  m=12  nonuniform  intervals  for  the  p 
coordinate.  The  consecutive  interval  length  ratio  was  q=1.2  starting  from  the  hole  edge.  The 
“fine”  subdivision  contained  ns=4  sublayers  per  ply  thickness  and  p  coordinate:  m=24, 
q=1.4. 

Asymptotic  analysis  revealed  that  a  real,  near-singular,  root  1  <  A,<1.15  is  present  at 
the  filled-hole  edge  in  this  problem.  The  singular  (0<A,<1)  and  the  near  singular,  root 
1  <A,<1.15  are  shown  in  Figure  12.  The  singular  root  of  the  open-hole  solution  is  also  shown 
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Figure  12.  Roots  of  the  Asymptotic  Soiution  for  Stress  Distributions  near  the  Fiiled-  and 
Open-Hole  Edge  Singularity  for  [45/-45]  Interface. 

for  comparison.  It  must  be  noted  that  X.=l  is  also  a  root;  however,  the  stress  field  which  it 
creates  is  constant  and  therefore  its  spline  approximation  projection  will  be  exactly  equal  to 
it,  i.e.  -  5”  =  0 .  It  was  found  that  at  locations  0=45, 135,  225,  315,  the  singular  root 

V  V 

becomes  equal  to  unity  as  well  as  the  near-singular  root.  At  these  locations  the  singularity 
disappears,  and,  according  to  the  fact  that  -s^-=0,m  =  1,2,3 ,  an  instability  of  determination 

of  Km  is  expected.  Two  roots,  A,i  <  1  and  1  <  X3<1. 15,  were  included  into  approximation  (15) 
at  each  circumferential  location  93<0<267  which  corresponds  to  the  contact  zone. 

Collocation  points  were  po  and  pi.  The  multiplicative  term  Ki  as  a  function  of  0  obtained 
with  the  two  subdivisions  is  shown  in  Figure  13.  The  dashed  regions  surround  the  0=135  and 
225  locations  where  the  root  equality  occurs.  Outside  these  regions  the  results  obtained  with 
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Figure  13.  Multiplicative  Factor  of  the  Singular  Term  of  the  Filied-Hole  Asymptotic 
Solution  in  the  Bearing  Loading  Problem  for  [45/-45]g  Laminate. 

the  two  subdivisions  are  very  close  and  indicate  convergence.  The  interlaminar  shear  stress 
component  is  shown  in  Figure  14  for  9=180  on  the  interface  \|/=0.  Dashed  lines  show  the 
stresses  (16)  at  the  interface  calculated  from  displacements  (7)  directly  using  Hooke’s  law. 
The  stress  values  in  the  upper  and  lower  ply  are  discontinuous  for  both  subdivisions.  The 
region  of  discontinuity  shrinks  while  increasing  the  subdivision,  however,  the  amplitude  of 
the  discontinuity  at  the  hole  edge  increases  at  the  same  time.  A  similar  effect  at  the  open- 
hole  edge  was  explained  due  to  directional  nonuniqueness  of  the  asymptotic  behavior  [2]. 
The  hybrid  approximation  (15)  provides  traction  continuity  and  converged  stresses  even  for 
the  coarse  subdivision,  as  shown  by  solid  lines.  Chaotic  behavior  of  the  multiplicative  factor 
Ki  near  6=135  and  225  has  been  anticipated.  However,  the  stress  behavior  at  these  locations 
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Figure  1 4.  interlaminar  T ransverse  Shear  Stress  in  the  0=1 80°  Cross  Section  in  the  Bearing 
Loading  Problem  for  [45/-45]s  Laminate. 


is  expected  to  exhibit  convergence  since  the  singularity  disappears.  Figure  15  displays  the 
same  as  the  previous  figure  at  0=135.  Only  the  stress  values  calculated  from  displacements 
(7)  directly  are  shown.  The  traction  discontinuity  observed  for  the  coarse  subdivision  is  an 
order  of  magnitude  lower  than  at  0=180.  It  cm  be  also  seen  that  the  fine  subdivision 
provides  traction  continuity  at  the  hole  edge,  thus  manifesting  absence  of  the  singularity.  It 
should  be  noted  that  further  development  is  required  to  describe  the  behavior  of  the 
multiplicative  factor  Ki  in  the  regions  where  the  singularity  impairs. 
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10.  CONCLUSIONS 


1 .  The  method  of  superposition  of  hybrid  and  displacement  approximations  was  extended  to 
provide  accurate  stress  fields  in  the  vicinity  of  the  ply  interface  and  the  hole  in  practical 
laminates  with  multiple  interfaces.  The  asymptotic  analysis  was  used  to  derive  the  hybrid 
stress  functions.  The  displacement  approximation  was  based  on  polynomial  B-spline 
functions. 

2.  The  coefficients  of  the  singular  terms  in  stress  solution  near  the  ply  interfaces  and  the 
open-hole  edge  were  determined  in  quasi-isotropic  [45/90/-45/0]s  laminates  under 
mechanical  and  thermal  loading.  Convergence  smdies  showed  that  accurate  values  of  the 
coefficients  of  the  singular  terms  can  be  obtained  with  the  coarse  out-of-  plane 
subdivision  of  one  sublayer  per  ply.  It  was  shown  that  for  laminates  with  multiple 
interfaces,  the  influence  of  singular  terms  on  adjacent  interfaces  is  important  for  coarse 
subdivision  convergence. 

3.  Coefficients  of  the  singular  term  of  the  asymptotic  expansion  were  determined  for 
[45/-45]s  laminate  under  bearing  loading  introduced  through  a  rigid  fastener. 
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APPENDIX 


All  =  Qf^cos^e  +  Q56^sin20  +  Q^g^sin^e,  A^  =  +  (Qg+Q^g^)sin0  cos0  + 

A22  =  Q66Cos20  +  Q^6^sin20  +  Q^§sin20,  A33  =  (^^^cos^Q  +  Q^5sin20  +  Q£^sm20, 

A2I  =Ai2  ,  A31  =Ai3=A32  =A23=0  , 

Bsi  =  (QS3V(^^5^)cose  +  ((^^3V(^5)sin0,  B32  =  (Qg+Qg)cose  +  (Qg+Qg)sin0, 

B2I  =Bi2=Bii  =B22=B33  =0 ,  Bi3  =B31,B23  =B32  , 

cii  =  (il  C12 = (H,  C22 = q2.  C33 = 

C2I  =Cl2  ,  C31  =Cll=C32  =C23=0  . 

Ply  stiffness  coefficients  Q^"'’  s=l,..,N  are  contracted  notations  of  the  fourth-order 
tensor  Cjf^ ,  used  in  the  text. 
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